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We introduce a new shell model of turbulence which 
exhibits improved properties in comparison to the stan- 
dard (and very popular) GOY model. The nonlinear 
coupling is chosen to minimize correlations between 
different shells. In particular the second order correla- 
tion function is diagonal in the shell index, the third 
order correlation exists only between three consecu- 
tive shells. Spurious oscillations in the scaling regime, 
which are an annoying feature of the GOY model, are 
eliminated by our choice of nonlinear coupling. We 
demonstrate that the model exhibits multi-scaling sim- 
ilarly to these GOY model. The scaling exponents are 
shown to be independent of the viscous mechanism as is 
expected for Navier-Stokes turbulence and other shell 
models. These properties of the new model make it 
optimal for further attempts to achieve understanding 
of multi-scaling in nonlinear dynamics. 



I. INTRODUCTION 

Shell models of turbulence Jj]-||] are simplified 
caricatures of the equations of fluid mechanics in 
wave-vector representation; typically they exhibit 
anomalous scaling even though their nonlinear in- 
teractions are local in wavenumber space. Their 
main advantage is that they can be studied via 
fast and accurate numerical simulations, in which 
the values of the scaling exponents can be deter- 
mined very precisely. Our interest in shell models 
stemmed from our efforts to develop analytic meth- 
ods for the calculation of the numerical values of 
the scaling exponents j(| . In trying to do so we dis- 
covered that the most popular shell model that was 
treated in the literature, the so-called GOY mod- 
els |l],|j , poses very tedious calculations because it 
exhibits slowly decaying correlations between ve- 
locity components with different wave-numbers. In 
addition it has large oscillations around the power 
law behavior in the scaling regime, making the nu- 
merical calculation of the scaling exponents less ob- 
vious than advertised. We therefore derived a new 
model which exhibits similar anomalies of the scal- 
ing exponents but much simpler correlation prop- 
erties, and much better scaling behavior in the in- 
ertial range. Since there is a significant number 
of researchers that are interested in this type of 



models independently of the analytic calculability 
of the exponents, we decided to present the model 
per se, discuss its good properties, display the re- 
sults of numerical simulations, and compare it to 
the standard GOY model. These are the aims of 
this paper. 

In Section 2 we review the popular GOY model, 
and explain the shortcomings that induced us to 
consider a new model. Section 3 introduces the 
new model, that we propose to call the Sabra 
model; we discuss the phase symmetries and cor- 
relations, stressing the much improved properties. 
Section 4 discusses numerical simulations from the 
algorithmic point of view. Section 5 contains the 
results of numerical simulations and fit procedures 
for accurate calculations of the scaling exponents. 
We believe that this section contains methods that 
should be used in the context of any shell model, 
and go beyond naive log-log plots. Section 6 
presents a discussion of the limitations in comput- 
ing high order exponents. We demonstrate that 
beyond Cs one needs exponentially long running 
times to extract reliable exponents. The evaluation 
of (io requires about one million turn-over times of 
the largest scales. We believe that similar limita- 
tions are important also in other examples of mul- 
tiscaling, including Navier-Stokes turbulence. Sec- 
tion 7 demonstrates the universality of the scaling 
exponents with respect to the viscous mechanism, 
and Section 8 offers a short summary. 



II. REVIEW OF THE GOY MODEL 

A. Basic Properties 

In the past considerable attention has been given 
to one particular version of shell models, the so- 
called GOY model . This model describes the 
dynamics of a complex "Fourier" component of a 
scalar velocity field that is denoted as u n . The 
associated wavenumber is 1-dimensional, denoted 
as k n . The index n is discrete, and is referred to as 
the "shell index" . The equations of motion read: 

du 

= i(ak n+1 u n+2 u n+1 + bk n u n+ iu n -i (1) 



1 



2 



+dc„_iU n _iM„_ 2 )* - vk 2 n u n + f n , 

where the star stands for complex conjugation. 
The wave numbers k n are chosen as a geometric 
progression 



k n = k X n 



(2) 



with A being the "shell spacing" parameter. /„ 
is a forcing term which is restricted to the first 
shells. The parameter v is the "viscosity". In the 
limit of zero viscosity one can arrange the model 
to have two quadratic invariants. Requiring that 
the energy 



(3) 



will be conserved leads to the following relation 
between the coefficients a, b and c: 







(4) 



A second quadratic quantity that is conserved is 
then 



H 



a/ c y 



(5) 



Our choice Eq. (g) is motivated by our reluctance 
to use the nonanalytic function \u n \. We will see 
that our definition yields £3 = 1 as an exact re- 
sult, similar to Kolmogorov's exact result for £3 in 
3-dimensional fluid turbulence. It was shown by 
numerical simulations that the choice of param- 
eters A = 2 and (a,b,c) = (1,-0.5,-0.5) leads 
to scaling exponents £ g that are numerically close 
to those measured in experimental hydrodynamic 
turbulence. 



B. Additional Properties 

The GOY model shares with Navier-Stokes tur- 
bulence an analog of the 4/5 law. Assuming sta- 
tionarity and using the quadratic invariants intro- 
duced above, we can obtain two identities involving 
third order correlations. Multiplying Eq. (|j) by u n 
we have, neglecting viscosity: 



^-S 2 (k n ) = 2k X n aXS 3 (k n+1 ) + bS 3 (k n ) (9) 
at 



+-^S 3 (k n -i) 



where 



Although non positive, this second invariant is 
often associated with "helicity" . 

The main attraction of this model is that it dis- 
plays multiscaling in the sense that moments of the 
velocity depend on k n as power laws with nontriv- 
ial exponents: 



(6) 



where the scaling exponents Q q exhibit non linear 
dependence on q. We expect such scaling laws to 
appear in the "inertial range" with shell index n 
much larger than the largest shell index that is ef- 
fected by the forcing, denoted as ti^, and much 
smaller than the shell indices affected by the vis- 
cosity, the smallest of which will be denoted as rid- 
We will refer to the moments as "structure func- 
tions" . For even q = 2m we use the usual defini- 
tion: 



S 2m (k n ) = (\u n \ 2m ) , 



(7) 



while for odd q = 2m + 1 we suggest the following 
definition: 

S2m+l(kn) = Im(M n _iU„U„ + i|u„| 2(m_1) ) , 

(GOY) . (8) 

The definition of the odd structure function differs 
from the usual definition S2m+i(k n ) — (\u n \ 2m+1 )- 



Pn 



2Re«/ n ) , 



(10) 



and obviously p n — for n > n^. In stationary 
conditions the rate of change of S^fcn) vanishes, 
and we find 

a\S 3 {k n+1 ) + bS 3 (k n ) + jS 3 (k n ^) = . (11) 
This equation has a solution in the inertial interval: 



S 3 (k n ) = — 

Kr>. 



A + B{- 

a 



(12) 



The unknown coefficients A and B can be found 
by its matching with the "boundary conditions" at 
small k n . To do so we can follow the considerations 
of Pissarenko et al Ji| and sum up Eq. (||) on all the 
shells from n — to an arbitrary shell M, where M 
is in the inertial interval. Using the conservation 
laws (i.e. a+b+c=0) we derive 



= 



2k M \aXS 3 (k M +i) + (b + a)S 3 (k M ) 



+ e. 



M 



2k 



fa\ M \ 

M^-J aXS 3 (k M +i) + (b + c)S 3 (k M ) 



(13) 

(14) 
+ S, 
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FIG. 1. Plots of the quantities z n = —k n S3(k„) with 
Sa(kn) taken as the stationary solutions (16). The 
fluxes e and 8 are related by e = cS/a and a = 1. 
The four panels have different values of c. Panel (a): 
c = 0.5, Panel (b): c = -0.5, Panel (c): c = 2, Panel 
(d): c=-2 



where the rate of dissipation e and the spurious 
mean 5 are defined as 



n=0 



TIL 

s>(!) 



(15) 



n=0 



Substituting the solution into Eqs. (jl|), (jjj) 
one relates the values of A and £? to those of the 
fluxes e and <5. Now Eq. (|r2[) becomes 



S 3 {k n ) 



1 



2fc„(a - c) 



5 



(16) 



There are four different types of functional de- 
pendence of S3(k n ) on k n , determined by the ratio 
c/a, as illustrated at Fig [lj For c/a < this func- 
tion has period-two oscillations which are caused 
by the existence of a non-zero flux of the second 
integral of motion which is not positively defined 
in this region. For c/a > the second integral is 
positively defined and the function is monotonic. 
For \c/a\ < 1 the role of the second flux becomes 
irrelevant in the limit n — > oo. Consequently the 
deviation of S^{k n ) from the scale invariant be- 
havior S 3 {k n ) oc l/k n decreases as n increases, see 
Panels (a) and (b) of Fig [l[ In contrast, in the 
case | c/a | > 1 the role of the energy flux becomes 
irrelevant in the limit of n — > oo. In this case the 
properties of the model are completely determined 
by the flux of the second integral, see Panels (c) 
and (d) of Fig [l| In the sequel we will focus on 
the region \c/a\ < 1. The reason for this is that 
Navier-Stokes turbulence never exhibits a region in 
which the energy integral is irrelevant. 

As we discussed, even in the "physical" region 
in which \c/a\ < 1 the sub- leading contributions 



(which are effected by the second integral) may in- 
fluence the apparent scaling behavior of the leading 
scale invariant contributions which are determined 
by the energy integral. In the region —1 < c/a < 
which is commonly discussed in literature, sub- 
leading contributions lead to period-two oscilla- 
tions which decrease upon the increase of n. These 
introduce additional problems in determining the 
scaling exponents. A simple way to eliminate this 
complication is to consider a "helicity-free" forc- 
ing chosen such that the flux of the second inte- 
gral ("helicity") would vanish identically. This is 
easilly achieved by selecting the forcing of the two 
first shells. From Eq. (|l5| ) we deduce 



6 = at cp o + api=0, P2=Pz 



-- . 
(17) 



For a random force which is Gaussian and 8- 
correlated in time 



{f n {t)f m it')) =al^ nm 5{t-t') 



one gets 



Pn 



(18) 



(19) 



For this type of forcing the condition of zero "he- 
licity" flux ([It]) is achieved by choosing the forcing 
to have the mean-square amplitudes 



c <7q + a o\ = 



(20) 



Under this condition the period-two oscillations 
disappear. 

The GOY model has some properties that make 
it undesirable for further analytic studies. It is 
best to exhibit these in comparison with the new 
(and we believe superior) model that we refer to 
as the "Sabra" model. 



III. SABRA MODEL: DEFINITION AND 
MAIN FEATURES 



A. Sabra model 

We propose the following equation of motion for 
the Sabra model: 

= i(ak n+ iu n+2 u* l+1 + 6fe n u n+1 u*_i (21) 
-cfc„_ 1 n n _ 1 'U n _ 2 ) - vk 2 n u n + f n , 

where for simplicity we assume that the coefficients 
a, 6, and c are real. As in the GOY model, conser- 
vation of energy in the inviscid limit is obtained if 
a+b+c= . 



4 



The fundamental difference with the GOY 
model lies in the number of complex conjugation 
operators used in the non linear terms. We show 
in the following that this slight change is responsi- 
ble for a difference in the phase symmetries of the 
two models. As a consequence, the Sabra model 
will exhibit shorter ranged correlations than the 
GOY model. Apart from this difference, all the 
calculations described in the previous section re- 
main valid. Both models share the same quadratic 
invariants and one can derive for the Sabra model 
another analog of the 4/5 law. We need to replace 
the definition of the odd orders correlators (|h ac- 
cording to: 

Ss(k n ) = Im(u„_iu„< +1 ) , (Sabra) 
S 2m +i(k n ) = lm{u n ^ lUn \u n \^ m -^u* n+1 ) . (22) 

Note that the shell index n is related to the inter- 
mediate shell involved in the correlation function. 



B. Phase symmetry and correlations 

Let us examine the phase transformation: 
u n ► u n exp(^2t/ n 

) • (23) 

The equations of motion of both the GOY and the 
Sabra models remain invariant under such trans- 
formations, provided that the phases 9 n are related 
by: 



_i=0, (GOY), 



T Cn T CrH 

9„_i + 6 n - 6 n+1 = , (Sabra) . 



(24) 



The phases 6 n can then be obtained iteratively 
from 9\ and 2l namely 

#i+3p = 0\, #3p+2 = #2, #3p = —8\ — 82 , (GOY); 

(25) 

8 n = -L[e 1 (ai- 2 -a n r 2 ) + e 2 (a 

v5 

a± = i(l =b \/5) , (Sabra) 



71 — 1 n—l\ 



(26) 



Although Eq. (|26|) has irrational numbers, it is easy 
to check that 



in = r n Vi + s„tf 2 , 



(27) 



where r„ and s n arc integer numbers which grow 
exponentially with n. 

Note that phases 9\ and 9 2 satisfy the equations 
of motion 



and they can be randomized by any small external 
forcing, ft means that any correlation functions 
which contain the phases #1, 9 2 or both phases 
must be zero. In our direct numerical simulations, 
see below, we confirmed that this is indeed the 
case. In the Sabra model there is only one nonzero 
2nd order structure function. All nondiagonal cor- 
relation functions vanish in the Sabra model 

S 2 (k n , fem) = (u n u* n ) = n/m (Sabra). (29) 

This is not the case for the GOY model for which 
there are correlations between shells separated by 
multiples of three: 



S 2 (k n ,k n+3p )^0 (GOY) 



(30) 



The relative simplicity of the Sabra model is seen 
also with regards to higher order structure func- 
tions. The Sabra model has only one non-zero 3rd 
order structure function S , 3(fc„) that couples three 
consecutive shells as defined by Eq. (p2|). All other 
3rd order structure functions vanish by averaging 
over the random phases 61 and 02- In contrast, in 
the GOY model there exists an infinite double set 
of nonvanishing correlation functions of 3rd order 
with given n. These are 



(u n u n+3p u n+ i +3q ) ± , (GOY) 



(31) 



The same phenomenon occurs also for higher order 
correlation functions. In the Sabra model the num- 
ber of non-zero correlation functions with finite n 
is much smaller then the corresponding functions 
in the GOY model, making it more convenient for 
theoretical analysis. 

To conclude this section we formulate a "con- 
servation law" that determines which correlation 
functions of the Sabra model are non-zero. Intro- 
duce a quasi-momentum K n for n-shell by 



where a is the golden mean, a 2 



(28) 



(32) 

li + 1. One can 
check that in the Sabra model the only non-zero 
correlation functions satisfy the following conser- 
vation law: the sum of incoming quasi-momenta 
(associated with u) is equal to the sum of outgoing 
quasi-momenta (associated with u* ). 



C. Additional properties 

In this subsection we show that the Sabra model 
exhibits the properties of the GOY model which 
were revealed in Subsect. p H . 

With this aim we compute from Eq. ( ^l| ) the 
time derivative of S^fcn, t): 



5 



dt 

= -21m 



dt 



(33) 



ak n {u* n u* n+1 u n+2 ) + bk n {Un_ x UnU 



n+1) 



-cfc n _i(w„_ 2 Wn-lW* 



2vk n (u n u n ) +p n 



where the forcing contribution p n was defined in 
Eq. 0. 

With the definition ( p2| ) of <S3(fc n ) this translates 
to the balance equation (^) derived for the GOY 
model. Note that these two models differ in the 
definitions of S 3 (k n ): Eq. (§) for the GOY model 



and Eq. (22) for the Sabra model. Clearly, S^AVi 



in the Sabra model has the same form ( J16| ) as in 
the GOY model and all the featu res of the GOY 
model discuused in Subsect. II B are relevant for 
the Sabra model as well. In particular one may 
eliminate the period-two oscillations by a proper 
choice (0) or ( po|) of the forcing. 

The reader should note however that in the case 
of the GOY model the second and the third order 
structure functions have additional long range cor- 
relations which do not appear in the balance equa- 
tion. This is a sickness of the GOY model that 
is eliminated in the context of the Sabra model, 
where what you see is what exists. Note also that 
the long range correlations for the GOY model ex- 
ists between shells separated by multiples of three 
(see, for example Eqs. (|30|), (|3l|)). These correla- 
tions are responsible for period-three oscillations 
in scaling plots of the GOY model. These annoy- 
ing oscillations are absent in the Sabra model by 
construction. Thus after elimination of the period- 
two oscillations (using "helicity-free" forcing) one 
finds scale invariant behavior of the structure func- 
tions almost from the very beginning of the inertial 
interval. 



IV. ASPECTS OF THE NUMERICAL 
INTEGRATION: STIFFNESS, FORCING 
AND DISSIPATION 

The numerical investigation of the Sabra model, 
as of any other stiff set of differential equations, 
calls for some care. We therefore dedicate this sec- 
tion to a discussion of the issues involved. A reader 
who wishes to consider the results only can skip 
this section and read the next one. 



A. Stiffness 

The main difficulty in integrating a shell model 
stems obviously from the stiffness of the system 
i.e. we are concerned with a wide range of time 



scales in the system. Within the inertial range, 
the equation is dominated by the non linear terms 
so that the natural time scale (in the Kolmogorov 
approximation) of the nth shell scales as: 



f 



f 



,2/3 



(34) 



Within the viscous range however, the dominant 
term is the viscous one and if the nth shell lies in 
this subrange, its natural time becomes: 



I 

vk! 



(35) 



We can now estimate the global stiffness of the 
system by quoting the ratio of the extremal time- 
scales: 



n 

T N 



TN 



2/3 



kd 



(36) 



x X 2[N+2(N-n d )-l]/3 



The global stiffness of the system thus depends 
both on the total number of shells N and on the 
width of the viscous region. Most of the results 
published in the literature are obtained with 22 
shells, a forcing restricted to the first shell and a 
viscous boundary beginning about the 18th shell. 
In this typical case, we have ti/tn ~ 6.6 x I0 5 . In 
this paper we typically use N — 34 with about 
6 shells in the viscous range. For this choice 

Ti /t n ~ 10 9 . 

To deal with this stiffness we chose from the 
library SLATEC the backward differentiation 
routine DDEBDF This routine is specially 

dedicated to very stiff problems. Although rather 
fast, its precision is not exceptional and it is rather 
sensitive to functions which are not sufficiently 
smooth. In cases of failure of the backward differ- 
entiation routine, the code switches automatically 
to a 4/5th order Runge-Kutta algorithm. Both 
routines adapt their step-size to fulfill a prescribed 
precision requirement. The backward differentia- 
tion routine adapts in addition its order between 1 
to 5. 



B. Random forcing 

We generate the random forcing to guarantee 
zero mean value of the velocity. We use a time 
correlated noise, with a correlation time chosen to 
be the natural time scale at the forcing shell: t = 



l/(k n 



Denoting the forcing term /, in case 



of an exponential correlation, the evolution of / is 
ruled by the equation 



G 



dt J 



h 77 . 

T 



(37) 



where 77 is an uncorrelated noise. The presence of 
this new equation in the system could in principle 
make the integration more cumbersome. Fortu- 
nately, the system being stiff, the typical time step 
used in the integration is very small compared with 
the forcing time scale r (six orders of magnitude 
in a typical calculation with 22 shells). This al- 
lows us to integrate f separately with a first order 
scheme. In the code, the forcing is updated at each 
new call of the integrator. The Gaussian exponen- 
tially correlated random forcing is computed (af- 
ter proper initialization) according to a first order 
scheme proposed by Fox et al j|]: 




f(t + dt) = f(t)E + cr v/-2(l- £ 2 )log(a) exp (i2w/3) 

(38) 

where E — exp(— dt/r), a is the standard devia- 
tion of / and a and (3 two random numbers be- 
tween and 1. 

C. Dimensional Analysis 

For the purpose of our numerical fits we con- 
sider, following jll|, the dissipative boundary rid, 
where the dissipative term balances the non linear 
term. At this boundary kdV% ld is of the order of 
vk\u. nd . In the viscous range n > rid one can guess 
a generalized exponential form: 



u n ~ k n exp 



where 



x = log; 



1 + V5 



(39) 



(40) 



We have studied the influence of the width of the 
viscous range on this exponential behavior. The 
results obtained for a system of 22 shells with var- 
ious viscosities are summarized in Fig. where 
we can see that the scaling behavior in the vis- 
cous range approaches slowly the asymptotic pre- 
diction. In the case of the largest viscosity used 
i/ = 8x 10~ 4 , we note that the asymptotic behav- 
ior starts at n ~ 15 while rid ~ 9. We can then 
consider that this width of 6 shells is the minimal 
one needed to properly describe the viscous range. 

In the inertial interval dimensional reasoning 
leads to K41 scaling: u n ~ (e/fc,,) 1 / 3 . This for- 
mula may be matched with (pE 



12.0 16.0 20.0 24.0 



FIG. 2. Modulus of u n in the Sabra model obtained 
by integration over 500 turn-over time scales with val- 
ues of the viscosities as shown in the figure. The dashed 
line represents the expected asymptotic behavior in the 
deep viscous regime. The slope of this line is given by 
equation (39). 



1+ _J_ 



exp 



(41) 

(/> 6 fc„J 1/8 - 
We will see that although the actual values of the 
exponents change due to multi-scaling, the form of 



where u nL ~ y/ f/k nL and kd 



the solution is rather close to reality, and Eq. (41) 
is a good starting point for numerical fits. 



V. NUMERICAL SIMULATIONS: RESULTS 

A careful determination of the scaling exponents 
is a delicate issue. With an infinite inertial range, 
we expect pure scaling laws. Despite its large size, 
the inertial range that we have in shell models re- 
mains finite. The most widely used method to de- 
termine the exponent is based on a linear regres- 
sion or on the determination of a local slope jl0| 
in log-log scale. In these methods one needs a cri- 
terion to choose the fitting range. The uncertainty 
in the scaling exponents comes obviously from the 
quality of the regression but also largely from the 
number of shells taken into consideration. We want 
to make the point here that these methods are not 
reliable, giving rise to a lot of confusion in the lit- 
erature. One needs to fit a whole function to the 
inertial and dissipative ranges simultaneously to 
achieve reliable estimates of the exponents in the 
inertial range. 

The definition of the scaling exponents can be 
matter of choice of the statistical object. Our pre- 



7 



(42) 




(43) 




9/3 ) 


sa 
o 



ferred definition is (fy and ( p2| ) for even and odd 
exponents respectiveiy. Two alternative choices 
were widely used in the literature, respectively 
based on the modulus of the velocity and on the 
energy flux: 



S q (k n ) = <|u»|«) , 
S q (k n ) = <|S„| 9/3 ) 

( |lm [a\u n u n+1 u* n+2 - cu n _iu n n* +1 



The latter definition allowed for a higher numer- 
ical precision in the context of the GOY model 
because the energy flux is not affected either by 
the genuine dynamical oscillation (due to to the 
hclicity flux) or by the period three oscillations. 
Beyond these different definitions of the statistical 
objects, we can also modify the definition of the 
scaling exponents themselves. In the framework 
of so called "Extended Self Similarity" , instead of 
writing S q (k n ) — Ak n ^ q one assumes a scaling re- 
lation between the structure functions of order q 
and of order 3: S q (k n ) = A[S 3 (k n )]^. 

These different definitions give a priori different 
sets of scaling exponents. An efficient compari- 
son is however difficult to set up in the case of the 
GOY model because of the various oscillations pol- 
luting the data. Moreover none of the techniques 
described so far took explicitly into account the 
finite size effects. The fitting procedure that we 
describe now is a first attempt to do so, and one 
of the results is that the exponents are universal, 
independent (for given parameters) of the choice 
of the statistical object. 

In light of the interpolation formula (|4l]), and en- 
couraged by the fact that the dissipative, stretched 
exponential behavior is rather nicely obeyed, we fit 
all our spectra to the following fit formula: 



F q (k n ) 



I 



k n 
* q ~k7 



>><i 



exp 



kn 
kd, q 



(44) 



This guarantees the right behavior at both asymp- 
totics. Note that we don't make any hypothesis on 
the form of the transition between the power law 
and the dissipative regimes. In fitting we minimize 
the following error function: 



£ = 



I - 



log F q (kn) 



log S q (k n ) 



(45) 



Here S q refers to the numerically obtained struc- 
ture function. We use the same fit formula for all 
the three definitions of statistical objects. The sum 



o.o 



-60.0 - 



-80.0 




in (45) was computed over the whole range except 



0.0 5.0 10.0 15.0 20.0 25.0 30 

log 2 (k n ) 

FIG. 3. Log-log plot of the structure functions 
S2(k n ) to S l 5(fc„) vs. k„ and of the results obtained 
using the fitting formula (44). The structure functions 
are represented by the symbols and the fits by the lines. 



the two first shells and the two last shells in or- 
der to limit the effect of the boundaries. It turns 
out that the minimum found in this procedure is 
sharp (as a function of £ q ) provided that we have 
a good fit of the S q over its whole range. To es- 
timate the relative error in the scaling exponents 
£ g we arbitrarily computed the values of C 9 that 
agree with values of £ that are twice the minimum 
value. These are the errors reported in all the ta- 
bles below. 

In all our simulations we used the parameter val- 
ues a = 1, b = c = —0.5 and ai/o- = 0.7. This 
choice eliminates the flux of helicity and corre- 
spondingly the period-two oscillations in the scal- 
ing plots. Typical fits for the structure func- 
tions from S2 to 55 for simulations with 34 shells 
(y = 4x I0 -11 , (To = 5x 10~ 3 ,) are shown in Fig. ||. 
In Table || we present the computed scaling expo- 
nents associated with three different definitions of 
g-order correlation functions. These results offer a 
very strong indication that the three scaling expo- 
nents of g-order correlation functions (with given 
q) are all the same. 

On the other hand, we can make the point that 
"Extended Self-Similarity" (ESS) |§ in its stan- 
dard usage does not seem be a useful approach in 
the present context for computing more accurate 
scaling exponents. In Fig. ^ we present S2{k n ) 
both as a function of k n and as a function of 
Sz{kn)- Even though superficially the ESS way of 
plotting seems to yield a longer linear plot, a care- 
ful examination shows a break in the inertial range 
scaling which occurs precisely at the crossover to 
dissipative behavior. We gain nothing from ESS 
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0.393 ± 0.006 


0.393 ± 0.007 


o 
Z 


U. r zU ± U.UU8 


U. / 2U ± U.UUo 


U. 1 iy ± U.UU / 


3 


1.000 ± 0.005 


1.003 ± 0.009 


1.000 ± 0.005 


4 


1.256 ±0.012 


1.256 ±0.012 


1.249 ± 0.003 


5 


1.479 ± 0.006 


1.488 ±0.013 


1.477 ± 0.004 


6 


1.706 ±0.015 


1.706 ±0.015 


1.691 ± 0.006 


7 


1.901 ±0.010 


1.910 ±0.020 


1.893 ±0.010 



TABLE I. Summary of the scaling exponents com- 
puted with a model of 34 shells 



-70.0 1 1 1 1 1 1 1 1 1 1 1 1 1 

-10.0 0.0 10.0 20.0 30.0 40.0 50.(1 

log 2 (k n ) 

FIG. 4. Log-log plot of S2{k n ) vs. k n and vs. 5*3 (fc n ) 
respectively. This plot shows that at least for this 
model, and when the accuracy is sufficiently high, 
ESS is quite useless in increasing the effective range 
of power law behavior. 

for this model. 

Nevertheless, for the limited aim of computing a 
precise value of £2, we can make use of the ESS idea 
provided that we fit the whole range. To do this 
we have to impose additional information on the 
fitting function. For the case of £2 we can employ 
the information contained in the balance equation 
(||), closing it with the ansatz 

S 2 (k n ) = A 2 \S 3 (k n )\^ . (46) 

Using ( ^6| ) and introducing z n = —k n S 3 {k n ), we 
can rewrite (^|) as: 

/ \2-C2 

Zn+l = Zn-1 +Hz n -l - Zn) - (/c„/fcj Pnl , 

(47) 

where K = (ja4 2 ) 1/(C2-2) and a = 1. 

Given zq and z\ , one can iteratively calculate z n 
and, consequently, S 3 (kn) and S 2 (k n ) in the range 
of k n , for which the ESS ansatz is valid with rea- 
sonable accuracy. Assuming for simplicity Zq = Z\ , 
the values of z n arc defined by 3 free parameters: 

As an example, we applied this procedure to the 
numerical data calculated with a = 1, b = —0.5, 
and v — 4 x Id -11 . The values of fitting param- 
eters corresponding to the global minimum of the 
functional £ (||) are z = 0.00126, A 2 = 1.80, and 
(2 = 0.728. To estimate the accuracy of the chosen 
fit parameters we have studied the dependence of 
£ on the deviation 6(2, 3A 2 , and Szq from their 
optimal values with two other parameters fixed at 



the optimal values. As before we define the error 
bar for each parameter interval for which £ takes 
on values which are twice the value at the mini- 
mum. With this definition z = 0.00126 ±0.00002, 
A 2 = 1.80 ± 0.06, and (2 = 0.728 ± 0.002. 

The accuracy reached here is higher than in the 
procedures described above. Most of the errors 
in the fit appear from the crossover region from 
power law to exponential decay. The analytically 
calculated S 2 (k n ) and Ss(k n ) near the onset of the 
viscous range are very sensitive to the value of £2- 
Therefore employing an adequate fit in this region 
(which uses additional a priori information con- 
tained in the balance equation) allows one to be 
more accurate. Note that we do not have such 
simple balance equations for higher order correla- 
tion functions and therefore a generalization of the 
procedure for higher orders is not available. 

VI. TESTS OF THE STATISTICAL 
QUALITY OF THE NUMERICAL DATA 

In evaluating the scaling exponents £ q one has to 
make sure that the structure functions S q {k n ) are 
calculated properly. This means that (i) the aver- 
aging time is sufficient for the representative statis- 
tics, and (ii) the numerical procedure produces an 
accurate realization u n (t). 

A. The PDF test for the averaging time 

In intermittent statistics one may need to wait 
rather long times before the appearance of rare 
events which nevertheless contribute significantly 
to the statistics of g-order structure functions of 
n-shells. This issue was carefully discussed by Lev- 
eque and She [ pd| in numerical simulations of the 
GOY model. They considered the waiting time 
Tn, q which are needed to evaluate safely q order 
correlations of n shells. They argued that times of 
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the order of 5 x 10 9 turnover times of the n-shell 
are required for q « 15. 

In the beginning of this subsection we estimate 
analytically the waiting time T n q which is needed 
to observe, say, 100 events contributing to S2 q {k n ). 
This is done using the probability W n , q to observe 
one rare event in which the value of the velocity 
u n hits the range that contributes mostly to the 
statistics of S2 q (k n ). Denoting by r n the decorre- 
lation time on the nth shell we estimate 



T n , q ~ 100r„/W„, 9 



(48) 



The probability W n ,q may be related to the PDF 
of the velocity at the nth shell, P n (u). For the sake 
of this estimate we take P ra (it) as a stretched ex- 
ponential. We do not imply that this distribution 
function is realized in this model (in fact we know 
that it is not consistent with multiscaling) . We 
use it only for the sake of an order of magnitude 
analytical estimate of the waiting time. Consider 



P n (v) =Cexp[- 



(49) 



where v is dimensionless velocity v ~ u/uq, uq is 
a characteristic velocity, ~ S 2 (k n ) and C is a 
normalization constant. One computes S2 q (k n ) as 



S 2q (k n ) = u\ q \ v 2q P n (v)dv 



(50) 



The integrand in (J50j) has a maximum at v 
where 



(2q/S) 



1/8 



(51) 



From (|49|) we can estimate the probability that v 
will attain a value within an interval of order of 
yfq ~ 1 around v q , which was denoted as W n>q . 
This interval of v values contributes maximally to 
S2q- Namely, 



W n , q ~P n {v q )=Cexp[-2q/S\ 



Equation ( p2\ ) leads to the estimate 
T n , q ~ 100r n exp(2g/(5), 



(52) 
(53) 



where r n is a characteristic decorrelation time for 
n's shell. The time T n>q is exponentially large. For 
instance, for 8 = 1 and 2q = 10, the averaging time 
required for accurate measurement of Siq (k n ) is of 
the order of 
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In every figure results are presented for 6250 (solid line) 
and 625 (dashed line) turnover times Tz. Already Ss is 
not accurate even with the longer run. 



Tn 



100e 



2 x 10V n 



(54) 



Admittedly this evaluation is rather rough. 
More accurate evaluations should be based on 
the numerically computed probability distribution 
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FIG. 8. Same as figure 5 but for shell #16. The 
solid line represents a longer run of 2.5 x 10 6 Ti6, and 
the dashed line a shorter run of 2.5 x 10 s t\q. 



functions as done for the GOY model in |ll|]. We 
plot the numerical value of (u/uo) 2q P(u 2 /uq) ver- 
sus (u/u ) 2 and see how noisy is the region which 
gives the main contribution to S^. Such plots for 
the third shell are presented at Fig [| for two real- 
izations, one averaged over 625 and the other over 
6250 turnover times of this shell, T3. In panels a, 
b, c we show the integrands for S2, Si, and Sq. 
One sees that S% and Si can be evaluated reason- 
ably well even from the shorter run, while 56 can 
be computed only from the longer run. Panels d, 
e and f present the analysis for S% , S10, and S12 
correspondingly. The evaluation of Ss is question- 
able even when the long run is used; The results 
for 5*10, S12, etc. are meaningless even for the run 
of 6250 turnover-times. This run is too short for 
this purpose. The same analysis for shell #7 (in 
the bulk of the inertial interval), with two runs of 
4OOOT7 showing that the improvement of the long 
run is not sufficient (see Fig. ||). We can hardly 
compute Sg from the longer run. In the viscous 
end of the inertial interval (say for the shell #12) 
our runs was ten times longer (4 x 10 5 Ti2) and 
the results can be seen in Fig. [?]. Now Ss can be 
computed reasonably well, but S10 is still buried 
in noise. Higher order structure functions cannot 
be estimated at all. Lastly, in Fig. || we present 
results for the shell #16 which belongs to the be- 
ginning of the viscous subrange. Here we have even 
longer run of 2.5 x 10 6 Ti6, resulting in a marginal 
improvement in the ability to compute 5*10. 

For the evaluation of the scaling exponent ( q one 
needs to compute S q (k n ) throughout the inertial 
interval. It appears that we can determine scal- 
ing exponents up to & from runs whose duration 
is about 5000 (longest) turnover times. In order 
to find exponents up to Cs we need runs of mini- 
mal duration of 10 5 (longest) turnover-times. An 
accurate determination of the exponent £10 calls 
for runs of about one million turnover times! Note 
that this estimate is in agreement with the simple 
analytical formula ( |54| ) presented above. Note also 
that these conclusions may very well be applica- 
ble also for the analysis of experimental data of 
hydrodynamic turbulence. The scaling exponents 
with our choice of parameters in the Sabra model 
correspond to those of Navier-Stokes turbulence, 
and it is likely that the far end of the probability 
distribution functions is as hard to reproduce in 
experiments as in our simulations. Since very long 
runs are rarely available in experimental data, this 
should serve as a warning that stated numerical 
values of higher scaling exponents should be taken 
with great caution. 
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B. Test of the numerical procedure 

The averaging time is not the only factor affect- 
ing the quality of the numerical data. Since the 
time dependence of u n (t) is highly intermittent we 
need to test carefully the ability of the numerics to 
cope with this. We need to check that the statisti- 
cal characteristics of the process u n (t) obey the ex- 
act relations imposed on the correlation functions. 
A simple test can be built around the first equa- 
tion of the infinite hierarchy relating S q {k n ) and 
S q +i{k n ). Consider Eq. (||) relating 5*2 and S3. 
In the inertial range, where the viscous term may 
be neglected, the largest term on the LHS (pro- 
portional to c) is balanced by the two first terms 
on the LHS. In the viscous range, where Ss(fc„) 
drops to zero very quickly, this term is balanced 
by the viscous term on the RHS. It is thus useful 
to rewrite Eq. (||) in the form of a "balance coef- 
ficient " (keeping in mind that S3 (fc n ) is negative 
and S 2 (k n ) is positive): 



v 



2) 



a \^3 {k n +i)\k n +i + b\S 3 (k n )\k n - vk 2 n S 2 {k n 



c|S , 3 (fc„_i)|fc„_i 



(55) 



If the numerical data satisfies the balance equation 
(||) accurately, the coefficients Cr? has to be unity 
for all n. In Fig. ^| we show that in our simula- 
tions this relation between S^^kn) and S 2 (k n ) is 
obeyed with accuracy better than 0.1%. However, 
this does not mean that less frequent events which 
contribute to higher order correlation functions are 
also correctly reproduced. To check the statistical 
reliability of S4(fc„) one can use the second equa- 
tion from the hierarchy, which connects S4 (k n ) and 
S^fcn) and so on. To measure this accuracy one 

(2") 

can define, analogously to Ci , a generalized bal- 



ance coefficient C, 



(2?) 



To define it we consider the 



time derivative of S 2q (k n )' 



dS 2q (k n ) 
dt 



= — 2qhn 



ak n ^ 



+bk n (^u* n _ 1 u* l u n+ i\u n \ 2< - q ^ 
-cfc„_i (u n - 2 Un-iu* n \u n \ 2{ - q ~ 1) 



T,+ \U n+ 2\ 



Iqvkl 



|2(9-1) 



(56) 



|2«A 



In the stationary case this gives: 

aS 2q+1 (k n +i)k n +i + bS2 q +i{k n )k n 

~^~ c ^2q+l{kn—l)k n —i = vk n S 2 q {k n ). 



(57) 



Here S 2q + 1 is defined by Eq.(g^) and we have in- 
troduced two additional structure functions: 

St q +i(k n ) = lm(u n -iu n u* n+1 \u n ± 1 \ 2< - c '^ 1 ^ . (58) 
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FIG. 9. Balance coefficient C. 



(2) 



for 22 shells, 



-0.5, average over 2500 largest turnover times. 



tj , in the form 



One can rewrite ( |57| ) , similarly to 
of a balance coefficient: 

Cl 29) (59) 
_ a S2q+i(k n +i)k n+ i + bS 2q +i{k n )k n + vk1S 2q {k n ) 
c \S 2q -\-i{k n -i)\k n -i 

Again, if the numerical data reproduce the balance 
equation ( |57| ) , the coefficient Cn^ has to be unity 
for all n. Testing this fact should be an integral 
part of the numerical solution of this model and 
similar models in the future. 



VII. UNIVERSALITY WITH RESPECT TO 
HYPERVISCOSITY 

"Hyperviscosity" in shells models amounts to 
changing the viscous term in ( ^l| ) with a term vk 2m 
with m > 1. The effect of hyperviscosity on shell 
models is a matter of controversy. It has origi- 
nally been argued by She and Leveque |l(| that 
in the GOY model there was no universality of 
the scaling exponents, the value of the latter be- 
ing strongly dependent on the dissipation mech- 
anism. The same observation has been made by 
Schorghofer |l3| et al and by Ditlevsen [Q . If true 
this observation would cast a doubt either on the 
relevance of shell models in turbulence studies or 
on one of the most widely accepted hypotheses in 
fluid turbulence: the universality of the exponents 
in the scaling range. Note for example that many 
direct simulation of 3D turbulence use hyperviscos- 
ity. On the other hand, Benzi et al have showed 
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FIG. 10. Coefficients Ci 2) [see Eq. (55) ] obtained 
for a model of 34 shells integrated over 250 forcing 
turnover time scales with dissipative terms propor- 
tional to k^, fc* and fc® respectively. 



Jig ] that in the case of shell models with eddy vis- 
cosity, the inertial exponents were independent on 
the particular definition used for the eddy viscos- 
ity. We made recently the point |D| that within 
the GOY model this phenomenon is nothing but 
a finite size effect which disappears when one in- 
creases the size of the inertial range. We dedicate 
this Section to showing that the same is true for 
the Sabra model. 

Before discussing the results we need to test 
our simulations for accuracy of the evaluation of 
the structure functions. To this aim we present 
in Fig. [Toj the balance coefficient C'n ^ (cf. Sect 



VI B ) for m = 1, 2,3. We tested the accuracy in a 
relatively short run of 250 forcing turnover times 
scales. The results indicate that even for this short 
run the accuracy of determination of the two low- 
est order structure functions is about 0.1% in the 
inertial range, but only about 1% in the dissipative 
range. Note that hyperviscosity makes the deter- 
mination of the structure functions in the viscous 
range (starting with the crossover region) some- 
what less accurate. In order to reduce the source 
of uncertainty and without loss of generality, we 
measured the exponents from the flux based struc- 
ture functions fl43|). 

In the following, we focus on the second and 
third order structure functions, using runs of dura- 
tion 1500 forcing turnover times scales, and offer a 
careful calculation of their apparent scaling expo- 
nents as a function of the number of shells used in 
the simulations, and of the order of the hypervis- 
cosity term m. We will show that the hyperviscous 
correction affects a finite number of shells in the 
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FIG. 11. Log-log plots of fc n S3(fcn) vs. k n in case 
of hyperviscosity of index m — 2 with different num- 
bers of shells and viscosities. The collapse has been 
obtained by shifting the abscissa. The solid line shows 
the constant behavior expected theoretically. One ob- 
serves clearly that the departure from this constant 
value only occurs in a region of about ten shells near 
to the viscous transition. When the inertial range is 
large enough, the predicted behavior is recovered. 



vicinity of the viscous transition. This number is 
relatively large, about ten shells or three decades of 
"length-scales" . The reason for this large effect is 
that we have a discrete model in which each shell 
interacts with 4 nearest neighbors. This means 
that with the standard shell spacing parameter 
A = 2, the local interactions spread over more than 
one decade of length scales. Nevertheless, we show 
now that this number remains unchanged when we 
increase the size of the inertial range, indicating a 
mere finite size effect. 

To see this point examine Figs, [ll] and |l2| 
in which we superpose results for fc ra S3(fc„) with 
m = 2 and m — 3 respectively, which were ob- 
tained in eight different simulations as detailed in 
the figures. The plots are as a function of log(fc„) 
with an appropriate shift in the abscissa. We see 
that in all cases the region of deviation from a con- 
stant function, associated with the theoretical ex- 
pectation Eq. (g) is of constant magnitude and of 
constant extent, independent of v or the total num- 
ber of shells. This is a clear indication that when 
the number of shells increases to infinity the scal- 
ing exponent £3 = 1 will be observed in a universal 
manner. 

Another way to reach the same conclusion is ob- 
tained by fitting structure functions as explained 
in section 6 to the formula (fl4|). We ran simula- 
tions for vn — 1,2 and 3 with N = 22, 28 and 34. 
The exponent x of the viscous tail for m — 2,3 ex- 
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FIG. 12. Same figure as above in the case of hyper- 
viscosity of index m = 3. Note that the amplitude of 
the bump is larger than in the previous case. 



hibits significant departures from its dimensional 
expectation (|40|). We obtained x ~ 0.75 for m = 2 
and x ~ 0.90 for m — 3, while a; ~ 0.69 for m = 1. 
These values which seem to be independent on the 
order of the structure function have then been used 
in the fitting procedure. The results for £2 and £3 
with normal viscosity were quite independent of 
N. On the contrary, hyperviscosity caused an ap- 
parent change in scaling exponents. However, as 
can be seen in Figs. O and 
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these values can be 
plotted as a function of l/[log(fc,i/fci)] 2 and they 
converge, for kd — > 00 to the values obtained for 
rn = 1. Note that log(kd/kx) is precisely the length 
of the incrtial interval, and kd was obtained from 
the fit. 



VIII. SUMMARY 

We presented a new shell model of turbulence, 
and demonstrated its improved properties in terms 
of simpler, shorter range correlations. The model 
exhibits anomalous scaling similarly to the GOY 
model and to Navicr-Stokes turbulence. In the fu- 
ture we will argue that the improved properties 
of this model help considerable in seeking analytic 
methods for the calculations of the scaling expo- 
nents. We used the opportunity of the introduc- 
tion of this model to examine carefully issues like 
the accuracy of determination of scaling exponents 
and the minimal length of running time required 
to achieve accurate structure functions. These con- 
siderations are model independent and pertinent to 
other examples of multiscaling as well. Lastly, we 
demonstrated the universality of the scaling expo- 
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FIG. 13. The apparent scaling exponent £3 as a 
function of the square of the inverse extent of the in- 
ertial interval, for for m = 1 (circles), m = 2 (squares) 
and m = 3 (diamonds). The tendency towards ("3 = 1 
is evident. 
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FIG. 14. The apparent scaling exponent ("2 as a 
function of the square of the inverse extent of the iner- 
tial interval, for m = 1 (circles), m = 2 (squares) and 
m = 3 (diamonds). The tendency towards £2 as found 
for normal viscosity m = 1 is evident. 
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nents with respect to the type of viscous damping. 
This universality was questioned in the recent lit- 
erature but we showed here for the Sabra model 
and previously pl| for the GOY model that there 
is no reason to doubt it. 
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